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ABSTRACT 

We investigate the interaction between a giant planet and a viscous circumstellar disk by means of 
high-resolution, two-dimensional hydrodynamical simulations. We consider planet masses that range from 
1 to 3 Jupiter masses (Mj) and initial orbital eccentricities that range from to 0.4. We find that a planet 
can cause eccentricity growth in a disk region adjacent to the planet's orbit, even if the planet's orbit is 
circular. Disk-planet interactions lead to growth in a planet's orbital eccentricity. The orbital eccentricities 
of a 2Mj and a 3Mj planet increase from to 0.11 within about 3000 orbits. Over a similar time period, 
the orbital eccentricity of a IMj planet grows from to 0.02. For a case of a 1 Mj planet with an initial 
eccentricity of 0.01, the orbital eccentricity grows to 0.09 over 4000 orbits. Radial migration is directed 
inwards, but slows considerably as a planet's orbit becomes eccentric. If a planet's orbital eccentricity 
becomes sufficiently large, e > 0.2, migration can reverse and so be directed outwards. The accretion rate 
towards a planet depends on both the disk and the planet orbital eccentricity and is pulsed over the orbital 
period. Planet mass growth rates increase with planet orbital eccentricity. For e ~ 0.2 the mass growth 
rate of a planet increases by ~ 30% above the value for e = 0. For e > 0.1, most of the accretion within 
the planet's Roche lobe occurs when the planet is near the apocenter. Similar accretion modulation occurs 
for flow at the inner disk boundary which represents accretion toward the star. 

Subject headings: accretion, accretion disks — hydrodynamics — methods: numerical — planetary systems: forma- 
tion 



1. Introduction 

A striking property of the extrasolar planets known to 
date is the range of their orbital eccentricities, which is far 
wider than that of planets in the solar system.^ Eccentrici- 
ties are typically ~ 0.2-0.3, but very low and high values are 
also found (e.g., Marcy et al. 2005). A variety of explanations 
have been proposed to explain the eccentricities. The high 
eccentricity cases may owe their explanation to the presence 
of a distant binary companion (Wu & Murray 2003; Takeda 
& Rasio 2005). In fact, the large eccentricity (e = 0.93) of 
HD 80606b could be the result of a three-body interaction 
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process, known as "Kozai cycle" (Wu & Murray 2003). More- 
over, numerical experiments show that the observed distribu- 
tion of orbital eccentricities for e > 0.6 can be reproduced by 
assuming the action of a Kozai-type perturbation (Takeda & 
Rasio 2005). 

For more typical eccentricities, other processes are likely 
to be at work. Planet-planet interactions involving scatter- 
ing on dynamical timescales is a possibility (Rasio & Ford 
1996). However, numerical experiments indicate that interac- 
tions between equal mass planets would produce more isolated 
planets on low-eccentricity orbits than those observed (Ford 
et al. 2001). Secular interactions between planets is another 
possible means of eccentricity excitation (Juric & Tremaine 
2005). This mechanism assumes that the planets have ap- 
propriate initial separations for evolution to occur on secular 
timescales (~ 10^" years). 

Disk-planet interactions can also give rise to planetary ec- 
centricities (Goldreich & Tremaine 1980; Artymowicz 1992; 
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Goldreich & Sari 2003). The evolution of orbital eccentricity 
depends on a delicate balance between Lindblad and corota- 
tion resonances. If the planet is massive enough to clear a 
gap, Lindblad resonances promote eccentricity growth, while 
corotation resonances damp eccentricity. If the corotation res- 
onances operate at maximal efflciency, they dominate over 
Lindblad resonances by a slight margin, and the result is ec- 
centricity damping (Goldreich & Tremaine 1980). Two mech- 
anisms have been proposed to weaken the effects of corotation 
resonances and thereby provide eccentricity growth. The first 
mechanism relies on a large enough gap about a planet's or- 
bit to exclude the corotation resonances from the disk, while 
leaving a remaining Lindblad resonance, the 1:3 outer reso- 
nance. This mechanism requires a massive enough companion 
and/or low enough viscosity to open a wide enough gap. Such 
a mechanism has been demonstrated for binary stars (Arty- 
mowicz et al. 1991). The second mechanism relies on the 
delicate nature of corotational resonances in their ability to 
weaken (saturate) when the disk viscosity is sufficiently small 
and the planet eccentricity sufficiently large (Ogilvie & Lubow 
2003; Goldreich & Sari 2003). A small finite initial eccentric- 
ity, typically e ~ 0.01, is required for the eccentricity to grow 
by this mechanism. Moreover, a companion object on a cir- 
cular orbit can drive eccentricity in a circumstellar disk. This 
process is believed to occur in disks involving 10-20 Jupiter- 
mass (Mj) companions (Papaloizou et al. 2001) and binary 
stars (Lubow 1991a,b). The disk transfers eccentricity to the 
planet's orbit by exchanging energy and angular momentum. 
The disk-planet system undergoes a coupled eccentricity evo- 
lution. 

One goal of this paper is to determine whether the eccen- 
tricity growth by disk-planet interactions occurs. Simulations 
by Papaloizou et al. (2001) suggested that orbital eccentric- 
ity is excited when the mass of the embedded body is larger 
than about 10-20 Jupiter-masses, while lower mass compan- 
ions experienced eccentricity damping. By applying higher 
resolution and simulating over longer timescales, we aim to see 
whether eccentricity growth can occur for lower mass, plane- 
tary mass, companions. In addition, we are interested in the 
effects of eccentricity on planet migration. 

Accretion of gas onto a planet is likely affected by the 
planet's orbit eccentricity. Circumbinary disks surrounding 
eccentric orbit binaries undergo pulsed accretion on orbital 
timescales (Artymowicz & Lubow 1996). A similar process 
could occur with eccentric orbit planets. The mean accretion 
rate could also be modified by the orbital eccentricity. An- 
other goal of this paper is to investigate this accretion process. 

We apply high-resolution hydrodynamical simulations to 
investigate disk-planet interactions over several thousand or- 
bital periods. In Section 2 we describe the physical model. 
The numerical aspects of the calculations are detailed in Sec- 
tion 3. Results on the growth of the disk eccentricity for a 
fixed planet orbit are presented in Section 4. Similar results 
on the growth of disk eccentricity were recently reported by 
Kley & Dirksen (2006). We describe results on the planet's 
orbital evolution in Section 5. Results on the pulsed accretion 
are described in Section 6. Finally, our conclusions are given 
in Section 7 



2. Model Description 

2.1. Evolution Equations 

We assume that the disk matter and planet are in copla- 
nar orbits. In order to describe the dynamical interactions be- 
tween a circumstellar disk and a giant planet, we approximate 
the disk as being two-dimensional, by ignoring dynamical ef- 
fects in the direction perpendicular to the orbit plane (vertical 
direction). This approximation is justified by the fact that the 
disk thickness is smaller than the vertical extent of the plan- 
etary Roche lobes for the cases we consider. Comparisons 
between two- and three-dimensional models of Jupiter-mass 
planets embedded in a disk indicate that the two-dimensional 
geometry provides a sufficiently reliable description of the sys- 
tem (Kley et al. 2001; D'Angelo et al. 2003; Bate et al. 2003; 
D'Angelo et al. 2005). 

We employ a cylindrical coordinate frame {O; r, (f), z} in 
which the disk lies in the plane z = and the origin, O, is 
located on the primary. The reference frame rotates about the 
disk axis (i.e., the z-axis) with an angular velocity fl and an 
angular acceleration fi. The set of continuity and momentum 
equations, describing the evolution of the disk, is written in a 
conservative form (see, e.g, D'Angelo et al. 2002) and is solved 
for the surface density, E, and the velocity field of the ffuid, 
u. The turbulent viscosity forces in the disk are assumed to 
arise from a standard viscous stress tensor for a Newtonian 
fiuid with a constant kinematic viscosity, v, and a zero bulk 
viscosity (see Mihalas & Weibel Mihalas 1999, for details). 

A locally isothermal equation of state is adopted by requir- 
ing that the vertically integrated pressure is equal to p = E 
and that the sound speed, Cs, is equal to the disk aspect ratio, 
H/r, times the Keplerian velocity. In this study, we use a 
constant aspect ratio throughout the disk. 

The gravitational potential of the disk, "1>, is given by 



$ = 



GMp 



GMp 



(1) 



Tp are the mass and the vector position of the planet, re- 
spectively. The length e represents a softening length that 
is needed to avoid singularities in the gravitational potential 
of the planet. The third term on the right-hand side of equa- 
tion (1) accounts for the acceleration of the origin of the (non- 
inertial) coordinate frame caused by the planet. We ignore 
disk self-gravity. For the disks we consider, the Toomre-Q 
parameter never drops below 4 during the simulations. 

The orbit of the planet evolves under the gravitational 
forces exerted by the central star and the disk material. Non- 
inertial forces arising from the rotation of the reference frame 
also need to be taken into account. Therefore, the equation 
of motion of the planet reads 



Tp = - ^^^^* ^J^'^'^^ rp - n X (n X rp) - 2 n X T-p 



(2) 



where the angular velocity and acceleration vectors of the ro- 
tating frame are defined as fl = Slf and ft = Clz, respec- 
tively. The second, third, and fourth terms on the right-hand 
side of the equation are the centrifugal, Coriolis, and angular 
accelerations, respectively. 

The angular velocity of the reference frame relative to a 
fixed frame, fl = fl(t), is chosen so as to compensate for the 
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azimuthal motion of the planet. Details of how this is achieved 
can be found in D'Angelo et al. (2005). Since the azimuthal 
position of the planet is constant, a circular orbit reduces to 
a fixed point in the rotating frame. For a planet orbit with 
e > 0, the planet radially oscillates between the pericenter 
distance, a(l — e), and the apocenter distance, a(l + e). In this 
frame, the trajectory is a straight line of length 2ae with center 
at r = a. This scheme has the advantage that the planet 
does not drift across the grid in the azimuthal direction. It 
always maintains a symmetric azimuthal position with respect 
to the zone centers (where eqs. [1], [3], and [4] are evaluated) 
while it moves radially. This method helps reduce artificially 
unbalanced forces that act on the protoplanet (see discussion 
in Nelson & Benz 2003a,b). 

The last two terms in equation (2), which represent the 
forces per unit mass exerted by the disk material on the planet 
and the star, are 



(r - rp) AMu{r) 



A/d {\r 



.2^3/2 



and 



r dMu{r) 



(3) 



(4) 



These two vectors are included in equation (2) only when 
the disk back-reaction is accounted for and the protoplanet is 
allowed to adjust its trajectory in response to the disk torques. 
When these terms are applied, the integration in equations (3) 
and (4) is performed over the disk mass comprised in the 
simulated domain, Md- 

2.2. Physical Parameters 

In these calculations, the stellar mass, M*, is the unit of 
mass and the initial semi-major axis of the planet's orbit, 
ao, represents the unit of length. The unit of time, to, is 
defined such that 1/to — sj G (M, -f Afp)/a|]. When we refer 
to the "orbital period" or "orbit" as length of time, we actually 
mean the initial orbital period fo/(27r). To provide estimates 
of various quantities in physical units, we adopt A'/» = 1 A/q 
and ao = 5.2 AU, thus one orbit is ~ 11.9 years. 

In order to treat the strong perturbations induced by gi- 
ant protoplanets and limit the influence of the imposed ra- 
dial boundaries, we considered an extended disk whose radial 
borders are at rmin = 0.3 ao and rmax = 6.5 ao. The disk ex- 
tends over the entire 27r radians in azimuth around the star. 
In physical units, the disk models cover a region from 1.56 to 
33.8 AU. The mass of the disk, in the absence of the planet, is 
M-Q = 1.58 X 10~^ M, within the radial limits of the simulated 
region. We used a constant disk aspect ratio, Hjr = 0.05. 
The unperturbed initial surface density is axisymmetric and 
scales as r~^/^, which produces an unperturbed density, at 
r — ao, equal to 75.8 gcm~^. However, given the large plane- 
tary masses considered in this investigation, we also included 
an initial gap along the orbit of the planet that accounts for an 
approximate balance between viscous and tidal torques (e.g., 
Lin & Papaloizou 1986). The initial gap profile is based on 
equation (5) of Lubow & D'Angelo (2006). The initial gap 
width is modified by a factor of 1 -I- aoeo, in order to account 
for the planet's initial orbital eccentricity, eo. 

The initial velocity field in the disk is a Keplerian one that 
is centered on the star and corrected for the rotation of the 
frame of reference. In order to account for the effects due 
to turbulence in the disk, a constant kinematic viscosity, v, 



was used. In terms of the Shakura & Sunyaev parameter 
(Shakura & Sunyaev 1973), we have a = qo (ao/r)^''^, where 
ao = 4 X 10~^ (in physical units, v ~ 10^^ cm^ s~^). Although 
spatial variations and time fiuctuations consistent with the 
MHD turbulence are not included, this relation yields a mag- 
nitude of a that is in the range found in MHD simulations 
(Papaloizou & Nelson 2003; Winters et al. 2003; Nelson & 
Papaloizou 2004). The infiuence of viscosity was explored by 
performing a few calculations with other ao values (1.2 x 10""^ 
and 1.2 x 10"^). 

Three planetary masses were considered: 1 Afj,2A'/,j, and 
3 Afj (i.e., the mass ratio q = Mp/Mt, ranges from 1 x 10^'^ to 
3x10"^). The planets were set on initially circular or eccentric 
orbits about one solar mass stars. We examined configurations 
with initial eccentricities, eo, up to 0.4. A complete list is 
given in Table 1. 

At time t = G the planet starts from the pericenter posi- 
tion, while its azimuth, (^p, remains constant (in the rotating 
frame) and equal to tt throughout the calculation. In order to 
allow the disk to adjust to the presence of the planet, we im- 
pose two stages to the evolution. During the first phase, the 
planet's orbit is static and terms (3) and (4) are not included 
in equation (2). During the second phase, the protoplanet is 
"released" from the fixed orbit and is allowed to react to the 
disk torques via the full form of equation (2). In the mod- 
els presented here the first phase lasts until the release time, 
t — fris, which ranges from 1000 to 1200 orbits. The second 
phase lasts from several hundred to several thousand orbits. 

The value adopted for smoothing radius e (in eqs. [1] and 
[3]) resulted from numerical experiments in each orbital eccen- 
tricity configuration. The chosen value of e was the smallest 
that prevented the integration time-step of the hydrodynam- 
ics equations from getting shorter than ~ 10~^ orbits. In 
models involving 1 Mj and 2Afj planets, we set e — 0.1 -Rh, 
where Rh = Vp (Bailey 1972) is planet's Hill radius 

(or sometimes called Roche radius) . In models involving 3 Mj 
planets, we applied softening lengths between 0.12 7?h and 
0.2 _Rh. The latter value was used at the highest initial orbital 
eccentricities, eo = 0.3 and 0.4. We found that torques within 
the Roche lobe do not dominate the planet orbital evolution. 
Moreover, the smoothing radius does not significantly affect 
planet accretion. Thus, e does not likely play an important 
role in these calculations. 

For simulations that account for the disk torques on the 
planet, an additional approximation is made, which is de- 
scribed at the beginning of Section 5. 

3. Numerical Method 

The equations of motion of the disk are solved numerically 
by means of a finite-difference scheme that uses a directional 
operator splitting procedure. The method is second-order 
accurate in space and semi-second-order in time (Ziegler & 
Yorke 1997). Hydrodynamic variables are advected by means 
of a transport scheme that uses a piecewise linear reconstruc- 
tion of the variables with a monotonised slope limiter (van 
Leer 1977). High numerical resolution in an extended region 
around the planet is achieved by using a nested-grid tech- 
nique (see D'Angelo et al. 2002, 2003, for details) with fully 
nested subgrid patches, whereby each subgrid level increases 
the resolution by a factor 2 in each direction. Tests on the 
behavior of the nested-grid technique applied in a reference 
frame rotating at a variable rate Q = Q{t) are reported in 
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Table 1 
Initial orbital eccentricities. 



Mp eo 

(Mj) 0.01 0.1 0.2 0.3 0.4 

1 • • • • • • 

2 • 

3 • • • • • 



the Appendix of D'Angelo et al. (2005). The equation of mo- 
tion of the planet (eq. [2]) is solved by using a high-accuracy 
algorithm described in D'Angelo et al. (2005). 

3.1. Grid Resolution 

The excitation or damping of a protoplanet's orbital eccen- 
tricity depends on a delicate balance between Lindblad and 
corotation resonances (Ogilvie & Lubow 2003; Goldreich & 
Sari 2003, and references therein). To study this balance, it is 
necessary to resolve the width of all the resonances involved 
in the process. The locations of first order eccentric Lind- 
blad resonances reside in a region that ranges radially from 
approximately 0.6 a to 1.6 a. Furthermore, calculations on 
the saturation of isolated noncoorbital corotation resonances, 
performed by Masset & Ogilvie (2004), suggest that a min- 
imum resolution requirement to avoid spurious damping of 
eccentricity is that 

Ar/a< 4.1 ^C^keq, (5) 

where k is the azimuthal wavenumber and are coefficients 
of order unity given in Ogilvie & Lubow (2003). Although 
nested grids provide a linear resolution < 0.01 in a region 
~ 2.5 a X O.Stt, the efi'ects of these corotational resonances 
are not localized in azimuth, so the nested grids do not sub- 
stantially improve their overall resolution. Calculations that 
follow the orbital evolution of the planet were executed with 
grid systems GSl and GS2 (see Table 2 for a description of 
all grid systems employed in this study). According to equa- 
tion (5), the global radial resolution we apply Ar/a — 0.02 
could produce spurious damping in the outer disk for orbital 
eccentricities e < 3 x 10~^/(fcg), with fc > 1. In the 1 Af j 
case, simulations with 0.01 < e < 0.015 could undergo spuri- 
ous eccentricity damping for the most distant outer resonance 
k — 2 only. For e > 0.015, there is no spurious damping. In 
the 3Afj case, simulations with 0.003 < e < 0.005 could un- 
dergo spurious damping for k — 2 only. For e > 0.005, there 
is no spurious damping. 

3.2. Mass Accretion Procedure 

Accretion onto the protoplanet was simulated by remov- 
ing material within a distance of race ~ 0.3 -Rh. Mass 
is removed by means of a two-step procedure according to 
which the removal timescale, race, is 0.03 orbital periods for 
r — Tpl < racc/2 and 0.09 orbital periods in the outer part 
of the accretion region, racc/2 < \r — rp\ < race. Notice that 
Tace is about cqual to the Keplerian period around the proto- 
planet (i.e., in the circumplanetary disk) at \r — rp\ = race/2. 
Tests were carried out to evaluate the sensitivity of the pro- 
cedure to both parameters race and race . These tests indicate 



that, if race is reduced by a factor of 1.5, the average mass 
accreted during an orbit varies by less than 10%. Increas- 
ing the removal timescale by a factor of 5/3 only affects the 
accretion rate by 5% (see also Tanigawa & Watanabe 2002). 
We also checked whether the accretion parameters influence 
the orbital evolution of the planet. Using the same tests, we 
found no significant differences over a few hundred periods of 
evolution. 

We do not add the removed mass to the planet mass over 
the course of the simulations. Doing so would increase the 
planet mass by ~ 1 Afj for the simulations reported here. 

3.3. Boundary Conditions 

Boundary conditions at the radial inner boundary, r — 
Vmin, allow outfiow of material, i.e., the accretion fiow toward 
the central star. The inner boundary conditions exclude in- 
flow away from the star into the computational domain. Two 
types of boundary conditions were used at the outer bound- 
ary (r = rmax): reflective and non-reflective. With the first 
kind, neither inflow nor outflow of material is permitted at the 
outer border, which behaves as a rigid wall. Although rmax 
is much larger than the apocenter radii of the planets in the 
simulations, wave reflection was observed at the outer border 
in this case, especially when the planet orbit was eccentric. 

In order to lessen the amount of wave reflection at the 
outer disk edge, we also applied non-reflective boundary con- 
ditions, following the approach of Godon (1996, 1997). In 
these circumstances, boundary conditions are not directly im- 
posed on the primitive variables (i.e, E and u), but rather on 
the characteristic variables (i.e., the Riemann invariants of 
one- dimensional flows). The basic idea is to let outflowing 
(inflowing) characteristics propagate through the grid bound- 
ary. The correct propagation of the flow characteristics across 
the border depends on how accurately the adopted solution 
exterior to the grid approximates the exact solution. 

At both disk radial boundaries, the flow is assumed to be 
Keplerian around the central star. This choice could lead to 
some small-amplitude wave excitation at r = r,nax, since the 
fluid tends to orbit about the center-of-mass of the system, 
rather than around the central star (see also Nelson et al. 
2000). The effects of such waves were checked to be unim- 
portant as they tend to dissipate within a short distance from 
the outer disk boundary 

4. Disk Eccentricity 

4.1. Global Density Distribution 

The global surface density in the disk is plotted in Fig- 
ure 1 for different planet's masses and orbital eccentricities. 
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Fig. 1. — Surface density in a disk containing a 1 j (left) and a 3Afj (right) planet at t — fOOO orbits. The axes are in units 
of the planet's orbital semi-major axis, ao- The grey scale bars are expressed in units of 3.29 x 10^ gcm^^. From top to bottom, 
panels refer to the configurations with orbital eccentricities e = 0, 0.1, and 0.3. In each panel the planet (smaller circle) is at 
pericenter and located at {X, y) = (e — 1, 0). 
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Table 2 

Grid systems used in the simulations. 



Grid 


GSl 


GS2 


GS3 


GS4 


level 


Nr X 


Nr X 


Nr X N^ 


Nr X N^ 


1 


313 X 317 


313 X 317 


623 X 629 


623 X 629 


2 


264 X 264 


264 X 264 


524 X 524 




3 




404 X 404 







Note. — The grid system GSl achieves the highest res- 
olution (Ar/ao, A(j!>) = (0.01,0.01) in the region {r,(f)) G 
[0.4, 3.0] ao X O.SStt, which is azimuthally centered on the planet. 
The grid system GS3 resolve the same region with a reso- 
lution {Ar/ao,A(l>) = (0.005,0.005). The latter resolution 
is obtained with the grid system GS2 in the region (r, 0) £ 
[0.5, 2.5] ao X 0.657r (also centered on the planet in azimuth). 




Fig. 2. — Azimuthally averaged surface density in a disk containing a IMj (left) and 3Mj (right) planet for various values 
of orbital eccentricity: e = (solid line), e = 0.1 (short-dashed line), e — 0.3 (long-dashed line), and e — 0.4 (dotted-dashed 
line). The horizontal axis is in units of the planet's orbital semi-major axis ao. For the vertical axis, (E) = lO"** corresponds to 
32.9 gcm~'^. The time is about 1000 orbits and the planet is at pericenter. 
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Fig. 3. — Gap region around tlie orbit of a 3A'/j planet with fixed e = 0.4 at various orbital phases. The axes are in units of 
the planet's orbital semi-major axis oq. The grey scale bars are expressed in units of 3.29 x 10^ gcm~^. The planet's orbit in 
the inertial frame is represented by a dashed ellipse. The smaller solid circle indicates the position of the planet around the star 
(larger solid circle). The motion of the planet (and the panel sequence) is counter-clockwise. 



Left panels refer to configurations with Mp = IMj, and right 
panels refer to configurations with Mp = 3 Mj . The orbital 
eccentricity increases from top (e = 0) to bottom (e — 0.3). 
In all of the panels, the planet is at pericenter {rp/a = 1 — e). 
For sufficiently long evolutionary times, the inner disk (disk 
interior to the planet) is largely depleted because of the tidal 
gap produced by the planet and because the grid does not 
have sufficient dynamic range to cover regions close to the 
star where an inner disk would reside. The outer disk is tidally 
truncated at a radial distance that depends on both Mp and e. 
Figure 1 shows that the size of the truncation radius increases 
with planet mass and eccentricity, analogous to the case for 
circumbinary disks (Artymowicz & Lubow 1994). In addition, 
a wave or wake propagates in the outer disk. The disk trun- 



cation also be seen in Figure 2, which shows the azimuthally 
averaged surface density for various cases. The outer gap 
edge becomes less steep as the orbital eccentricity increases 
and does not change significantly over an orbital period, as 
illustrated in Figure 3. 

After a few hundred orbits, the disk region near the planet 
becomes eccentric, even if the planet is on a circular orbit. 
This effect can be seen in Figure 4, which illustrates the mo- 
tion (relative to the star) of the center-of-mass of the disk 
interior to r = 3 ao for different planetary masses and orbital 
eccentricities. 
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4.2. Mode Analysis 

We performed a mode decomposition of the surface density 
distribution of the disk adopting an approach along the lines 
of Lubow (1991b). We defined a mode component at each 
radial location in the disk as 

M\ijl - .T(E)a+M IS'' ^^''^ ^^'"^^^ '''^ ' 

where Sip = \J G (A/* + Mp)/a^, the angle 6 is the azimuth 
relative to an inertial reference frame, and 27r(E) — J^^'Edd . 
The time interval T is a ten-orbit period interval, beginning at 
a pericenter passage of the planet. The time integration is re- 
peated every interval T. The functions / and g are either sine 
or cosine. This decomposition corresponds to a Fourier trans- 
form in both azimuth and time. The amplitude (or strength) 
of the mode is 

[M'^:r'y+[K::rfY'- (7) 

In order to obtain the strength of a mode integrated over a 
radial interval [ri,r2], the mode components M^l'^-j are first 
averaged between ri and r2 and then substituted into Equa- 
tion (7). A tidally disturbed non-eccentric disk has modes 
present that all have k = I. For a disk ring to be eccentric, 
the mode strength associated with the pair {k,l) = (1,0), 
known as eccentric mode, must be non-zero, i.e., 5(1,0) > 0. 
To follow the eccentricity evolution, we analyze S(i_o) as a 
function of the time. 
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Fig. 4. — Motion of the center-of-mass of the disk region 
r < 3ao relative to the star (located at the origin). The disk 
center-of-mass starts at the origin, but is plotted at later times 
t > 300 orbits. The long-dashed and solid lines refer to the 
cases with Mp — IMj and Adp = 3Mj, respectively, with a 
fixed circular orbit planet. The short-dashed line is for the 
model with Mp — 3Mj and fixed planet eccentricity e — 0.1. 



4.3. Eccentricity Growth 

The amplitude of the eccentric mode, as defined in equa- 
tion (7), is shown in Figure 5 during the evolutionary phase in 
which the planet is kept on a fixed orbit. This Figure refers to 
the configurations with Mp — IMj, 2Mj, and 3Mj planets 
having e = 0. The disk region that undergoes a substan- 
tial eccentricity growth is that between r ~ ao and r ~ 2ao 
(Fig. 5, left panel). The relatively small initial value of ^(i.o) 
is likely due to a transient effect, as the initially circular disk 
adjusts to the presence of the planet. Farther away from the 
planet's orbit, the mode strength drastically decreases, and 
the disk is nearly circular (Fig. .5, middle and right panels). 
The eccentricity growth proceeds very rapidly during the first 
200 orbits and oscillates afterwards with some reinforcement 
in the 3 Afj case. 

The eccentricity driven in the region ao < r < 2ao by 
the 1 Afj planet (e = 0) is rather small compared to that 
driven by the other two planetary masses. However, with an 
initial orbital eccentricity e = 0.1, a IMn planet is able to 
sustain eccentricity growth in the disk, as indicated by the 
solid line in the left panel of Figure 6. In the 3 Afj case, an 
initial orbital eccentricity of up to e = 0.3 induces a larger 
amplitude eccentricity perturbation only at the beginning of 
the evolution. But the long-term evolution of ^(i^o) is not 
greatly affected (Fig. 6). 

The simulations for mode analyses were generally per- 
formed with the grid GS4. We also ran simulations and com- 
puted modes by using the grid system GSl and found results 
consistent to those obtained with the higher resolution grid. 

4-3.1. Influence of Viscosity 

The sensitivity of disk eccentricity growth to the disk kine- 
matic viscosity was examined for a system containing a 3 Mj 
planet on an initially circular orbit. Two additional values 
of ao (i.e., a at r = ao) were considered: ao = 1.2 x 10~^ 
and Qfo = 1-2 x 10~^, which are respectively a factor of 3 
smaller and larger than the standard value ao = 4 x 10^''. In 
the lowest viscosity model, the overall evolution of the ampli- 
tude of the eccentric mode closely resembles that of the model 
with standard viscosity (see Fig. 5, solid line), but was roughly 
20% larger. In the highest viscosity model, the eccentric mode 
strength reached a maximum of 0.3 at t — 200 orbits and then 
decays. Around t — 1000 orbits, ^(i^o) ~ 0.1 and continues 
to decline. The effects of viscosity were also investigated with 
a calculation involving a 1 Afj planet on a circular orbit and 
Qfo = 1.2 X 10~^. In this case, the mode amplitude S(i,o) is 
sustained at about 0.15 for t > 200 orbits, unlike the case 
of declining eccentricity at higher viscosity that is displayed 
as a short-dashed line in the left panel of Figure 5. These 
results suggest that disks with a less than a few times 10~^ 
experience sustained eccentricity, while disks with a in excess 
of « 10~ do not. The weakening of disk eccentricity with 
viscosity was also found in simulations by Kley & Dirksen 
(2006). 

4.4. Analytic Model 

In the case of superhump binaries, disk precession is dom- 
inated by the gravitational effects of the companion which 
causes prograde precession (Osaki 1985). On the other hand, 
pressure provides a retrograde contribution which is somewhat 
weaker (Lubow 1992; Goodchild & Ogilvie 2006). In the case 
of a circular orbit planet, the gravitational contribution to 
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Fig. 5. — Amplitudes of the eccentric mode, associated with the pair {k,l) — (1,0) (see eq. [7]), versus time for three different 
disk regions, as indicated in the top-left corner of each panel, for 1 Mj (short-dashed), 2Mj (long-dashed), and 3Mj (solid) 
planets. In all cases, the planet resides on a fixed circular orbit. The mode strength progressively weakens as regions farther 
from the planet's orbit are considered. 




Fig. 6. — Left. Mode strength ^(i.o) as function of time for models with Mp — 1 Mj of fixed planetary orbital eccentricity e = 
(dashed line) and e = 0.1 (solid line). Right. Sfi o) versus time for models with Afp = 3 Mj of fixed planetary orbital eccentricity 
e = (short-dashed line), e = 0.1 (long-dashed line), and e = 0.3 (solid line). 
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precession is expected to be weaker than the disk's pressure 
contribution. The magnitude of the pressure induced preces- 
sion rate is ~ (H/r)^ Q. For the disks simulated in this paper, 
the precession timescale is then ~ 10'^ orbits. We estimated 
the gravitational precession due to a 3 Mj planet for the case 
plotted in Figure 5 at 1000 orbits and find that this rate is 
about 10 times smaller than the pressure precession rate. We 
then expect the precession to be pressure-dominated and ret- 
rograde, with a timescale of about 10"^ orbits. This estimate 
in accord with the simulation results for a 3 Afj planet in a 
circular orbit in Figure 4. 

The group velocity for an eccentric mode is estimated by 
using the dispersion relation for an m = 1 disturbance in 
a Keplerian disk perturbed by pressure. The group velocity 
is Vg ^ {H/r)'^Q,r, where we assume the radial wavenum- 
ber \kr\ ~ 1/r and that the pattern speed is small compared 
to Q,. The group velocity leads to eccentricity propagation 
timescales of order (r/H)'^ ~ 10'^ orbits. This timescale is 
consistent with the localization of the eccentricity over course 
of the simulation of ~ 10** orbits to within a region of order 
2 a, as seen in Figure 6. Over longer timescales the eccentric- 
ity would spread further. 

We analyze growth of eccentricity in an outer disk that is 
perturbed by a planet on a circular orbit. Disk eccentricity 
growth via the 3:1 resonance for inner disks of superhump bi- 
naries occurs on a timescale ~ 0.1w/{rq^) binary orbit pe- 
riods, where w is the radial width of the eccentric region 
(Lubow 1991a). We now consider the case for outer disks 
perturbed by planets. For simplicity, we assume that the ec- 
centric corotational resonances in the disk are saturated (i.e., 
of zero strength). This situation is likely to hold if the disk 
eccentricity is of order 0.01, by analogy with the case of an 
eccentric orbit planet interacting with a circular disk (Ogilvie 
& Lubow 2003). Following the mode coupling analysis for ec- 
centric Lindblad resonances involving a circular orbit planet 
(Lubow 1991a), we find that the disk eccentricity growth rate 
associated with a particular eccentric outer Lindblad reso- 
nance is given by 
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where w denotes the radial extent of the eccentric region and 
prime denotes differentiation in r. All quantities are evaluated 
at the location of the eccentric Lindblad resonance associated 
with azimuthal wavenumber m, having Q = mflp/{m + 2) and 
r = Vm = [{m + 2)/m]'^^^ a. Quantities Um and Vm are the 
velocity components in a circular disk associated with poten- 
tial iprn, and 5 is the Kronecker delta function in the indirect 
potential term. 

The combined effects of all the resonances is determined 
by the local disk density at each resonance. The growth rate 
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Fig. 7. — Large dots: contributions to the growth rate sum 
for Equation (13) in units of Qp as a function of azimuthal 
wavenumber m associated with an eccentric Lindblad reso- 
nance. Small dots: normalized azimuthally averaged surface 
density at each eccentric Lindblad resonance, 1, as 

a function of m. The plot is for a 3 Mj planet with a density 
profile obtained from simulations at 500 orbits. 



of the mass- weighted eccentricity is then 
_ 2nw 

IVlc 



(13) 



where Mo is the mass of the eccentric region and the sum 
is taken over the active eccentric Lindblad resonances, and 
mmax — r/H due to torque cut-off effects (Goldreich & Tre- 
maine 1980). 

The contribution of each resonance in the above sum for A 
for a 3 Mj planet is shown in Figure 7. We consider times be- 
yond several hundred orbits, when the width of the eccentric 
region w ~ 2 a. Figure 7 shows that there is a weak contribu- 
tion from the outermost resonance, the 1:3 resonance, corre- 
sponding to m = 1. Even though the density at this resonance 
is the largest, a nearly complete cancellation occurs in Ai, due 
to effects of the indirect term in potential tpi. If only the 
m = 1 resonance were involved, then the eccentricity growth 
timescale would be very long, ~ 10"" orbits. The growth rate 
contributions from regions closer to the planet are weakened 
by the lower density, but strengthened by the larger values 
of Xm- At a time of 500 orbits, the density near the planet 
is small enough that the outermost resonances (2 < m < 5) 
provide most of the growth rate. At a time of 500 orbits, the 
eccentricity growth timescale 1/A is about 600 orbits, which is 
in very rough agreement with the average growth rate implied 
by Figure 5 for the innermost region, although there are con- 
siderable fluctuations in the simulations.'^ At earlier times. 



•^To compare the growth rate defined by A with simulations, it is 
better to adopt a similar mass-weighted eccentricity. This requires 
a slightly different definition of "Sji q) than given by Equation (6). 
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the growth rate is higher because w is smaller. We conclude 
that the disk eccentricity growth is possible over ~ 10^ orbits 
in the case of planets because of contributions from several 
resonances that lie in the disk edge/gap region. This situ- 
ation differs from the superhump binary case where only a 
single resonance is involved, since the disk extends relatively 
closer to the perturber in the planet case, due to its weaker 
tidal barrier. 

Increased viscosity affects the disk eccentricity in multi- 
ple ways. It leads to further disk penetration of the planet's 
tidal barrier, which could lead to stronger eccentricity growth. 
On the other hand, viscosity also acts to unsaturate (or 
strengthen) the corotation resonances which act to damp ec- 
centricity. Furthermore, the viscosity acts to damp the non- 
circular motions associated with the eccentricity. 

5. Planet Orbital Evolution 

At time t — t^is (between 1000-1200 orbits) the planet was 
allowed to adjust its orbit in response to the gravitational 
forces exerted by the surrounding disk material. We generally 
neglected torques on the planet due to gas within a distance 
of 0.5 i?H from the planet. However, analyses of the torque 
distribution at the release time indicate that torques within 
the Roche lobe do not dominate the orbital evolution of the 
planet. So we believe this procedure does not likely lead to 
major errors in the planet's orbital evolution. All calcula- 
tions presented in this section employed grid GS2, except for 
the 3Mj cases on initially circular orbits (eo = 0), which 
employed grid GSl. In order to analyze the orbital evolu- 
tion of the planet after release, we calculated the osculating 
elements of the orbit each few hydrodynamics time-steps (ap- 
proximately every 0.01 orbits). To remove short-period os- 
cillations, we computed the mean orbital elements (Beutler 
2005) by using an averaging period of one orbit. Throughout 
the paper, the planetary orbital eccentricities and semi-major 
axes in the simulations refer to the mean orbital elements. 

5.1. Orbital Eccentricity 

Simulations by Papaloizou et al. (2001) showed that the 
interaction between an initially circular disk and a circular 
orbit planet with mass > 20 Adj can lead to the growth of 
disk eccentricity and planetary orbital eccentricity. They also 
found that this interaction can be more efflcient at augmenting 
orbital eccentricity than direct wave excitation at the outer 1:3 
Lindblad resonance in a non-eccentric disk (e.g., Artymowicz 
1992) . We aim at determining whether a similar phenomenon 
can occur also in the Jupiter-mass range. 

Figure 8 shows that the interaction between a planet and a 
disk leads to orbital eccentricity growth for the 2 Mj (dashed 
line) and 3 M,j (solid line) cases. During the initial growth 
of e for the 3M,j planet, the rate is e ~ 1.3 x 10^** orbit^^. 
This value is ~ 1.6 times that exhibited by the 2Mj planet. 
The eccentricity growth stalls when e ~ 0.08 for both planet 
masses. The planets may be experiencing some variation in 
their eccentricity forcing due to the phasing of their eccen- 
tricities relative to the disks'. After 1500-1600 orbits from 
the release time, the orbital eccentricity starts to increase 



When we apply that definition, we obtain eccentric disk evolution 
in the 3Mj case that is similar to the 3Mj results in Figure .5, 
except that the magnitude of 5(i qj in the innermost region is about 
a factor of 8 smaller. The average growth rate of Sji g) in the 
innermost region is about the same over 1000 orbits. 
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Fig. 8. — Evolution of the (mean) orbital eccentricity of 2 Mj 
(dashed line) and 3 Mj (solid line) planets after the release 
time, tris = 1000 orbits. 



again with a growth timescale that is comparable to the ini- 
tial growth timescale, te = e/|e| « 2.3 x 10^ orbits, for both 
planetary masses. The average growth timescale is shorter 
than the standard Type II migration (or viscous diffusion) 
timescale (see § 5.2). Over the last 1000 orbits of the simula- 
tion, the eccentricity of the 3 Mj planet increases very slowly, 
at a rate e « 2 x 10~® orbit 

The model with Mp = 1 iV/j and initial zero-eccentricity 
(Fig. 9, dashed line) shows a much slower orbital eccentricity 
growth, reaching e = 0.02 after 3000 orbits from the release 
time. As described in Section 4.3, the eccentric perturbation 
induced by the planet on the disk is also rather weak compared 
to that excited by the 2 Mj planet. At the average growth rate 
e ~ 7 X 10~^ orbit"^, it would take on the order of the viscous 
diffusion timescale to reach e ~ 0.1. 

In order to evaluate to the extent of orbital eccentricity 
growth, we used configurations with fixed non-zero planet ec- 
centricities prior to release, eo. Figure 9 also shows the orbital 
eccentricity evolution of a Mp = IMj planet with eo = 0.01. 
After release, e oscillates about the initial value. The oscil- 
lation grows in amplitude and, during one of these cycles, e 
increases from to 0.09 within about 1300 orbits. In this 
case, Te is of order the viscous diffusion timescale. We sim- 
ulated several models with eo > 0.1 (not plotted) and found 
that there was generally a reduction of the orbital eccentric- 
ity, with some exceptions though. For example, in a model 
with a 1 Af,j planet and eo = 0.1, e underwent small ampli- 
tude oscillations about the initial value, with periods of a few 
hundred orbits. This occurrence may be related to the rela- 
tively large eccentricity driven in the outer disk. Some models 
with eo > 0.2 showed a rate of change of e that diminishes 
in time. In these cases the evolution was generally monitored 
for less than 1000 orbits. Longer time coverage simulations 
are required to assess the long-term behavior of these config- 
urations. 
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Fig. 9. — Orbital eccentricity versus time of Jupiter-mass 
models with different initial orbital eccentricities: eo = 
(dashed line) and eo = 0.01 (solid line). The release time 
is 1100 orbits. 
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Fig. 10. — Evolution of the semi-major axis of planets: Afp = 
1 M,j (short-dashed line), Mp — 2Mj (long-dashed line), and 
Mp = 3 Mj (solid line). The release time is about 1000 orbits. 
The initial eccentricity eo = in the 3 cases. The change in 
migration rates for the 2 Mj and 3 Mj cases at later times is 
related to their increased orbital eccentricity (see Figure 8). 

5.2. Radial Migration 

We describe here some results on the migration of eccentric 
orbit planets. We plan to explore this issue further in a future 
paper. Radial migration of planets in the mass range consid- 



ered in this study is expected to be in the standard Type II 
regime, which is characterized by an orbital decay timescale 
TM = ao/|d| = 2al/{3u) = (Ward 1997). The Type II 
migration rate depends only on the viscous timescale of the 
disk near the location of the planet and is independent of the 
disk density, provided the disk is locally more massive than 
the planet. Type II migration is based on the assumption that 
the gap, which separates the inner and outer disks, is devoid 
of material. In this case, the planet torques approximately 
balance the viscous torques at the gap edges. 

Figure 10 plots the evolution of the semi-major axes, af- 
ter the release time, of models with three planetary masses: 
Mp = IMj (short-dashed line), Mp = 2Mj (long-dashed 
line), and Mp — 3M] (solid line). The Figure shows that the 
initial migration rate depends on the planet's mass, which is 
inconsistent with the Type II prediction. We may expect some 
dependence of the migration rate on planet mass, because ra- 
tio of planet to disk mass is non-zero (Syer & Clarke 1995; 
Ivanov et al. 1999). In order to explore this result further, we 
have used a one-dimensional disk evolution code, along the 
lines of Lin & Papaloizou (1986). We used the torque density 
per unit mass given in equation (4) of Lubow & D'Angelo 
(2006). We checked that the results are insensitive to the de- 
tails of the torque density, provided that it is large enough to 
produce a gap. Increasing the torque density everywhere by a 
factor of 2 produced a small change in the migration rate (less 
than 1%). We adopted the same disk and planet parameters 
as in the two-dimensional simulations with zero planet eccen- 
tricity. In short, we find that the largest contributing factor 
to this non-Type II behavior is the lack of a substantial inner 
disk in the two-dimensional calculations. There is also some 
effect due to the non-zero planet-to-disk mass ratio. 

A comparison of orbital migration between one-dimen- 
sional and two-dimensional models is shown in Figure 11. 
In this comparison we used the azimuthal averaged surface 
density distributions in Figure 2 as initial conditions for the 
one- dimensional models. As seen in Figure 11, the one- 
dimensional (zero eccentricity) migration rates, for a very 
low density inner disk, agree well with the two-dimensional 
rates at early times after release, while the planet eccentric- 
ity is small. For undepleted initial inner disks, we find that 
one- dimensional models have about the same migration rates 
for these two planet masses. It is possible that the two- 
dimensional simulations have an inner boundary Vmin that 
is too large to resolve the inner disk. More complete zone 
coverage of the inner region in the two-dimensional calcula- 
tions might reveal an inner disk that acts to make the mi- 
gration rate less dependent on mass, as indicated by the one- 
dimensional simulations. In spite of these possible limitations 
of our two-dimensional simulations, we describe below some 
interesting aspects of the migration of eccentric orbit planets 
in two-dimensional disks. 

As a planet's orbital eccentricity grows toward values of 
about 0.08, the rate of migration slows significantly (see 
Fig. 10). Over the last 1000 orbital periods of the calcu- 
lated evolution, the 3 Mj planet exhibits a migration speed 
a ~ —2 X 10~® ao per orbit, with a tendency towards further 
reduction. This migration rate is about a factor of 30 smaller 
than the rate at release time. The 2 Mj shows an even more 
drastic reduction of the migration rate that actually reverses 
and becomes positive around t — tris « 4200 orbits. The out- 
ward migration speed is d « 1 x 10~^ ao per orbit at the end 
of the simulation. The migration speed of the 1 Mj planet 
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Fig. 11. — Orbital migration of a 1 Mj and a 3Mj planet 
according to one-dimensional simulations (dots) and two- 
dimensional simulations (lines) that use the same disk and 
planet parameters. The two upper curves are for the 
3M,j planet. The one-dimensional simulations have an ini- 
tially depleted inner disk, whose density distribution matches 
that of the two-dimensional simulations at planet release. 
Note the good agreement between one-dimensional and two- 
dimensional simulations for the 1 Afj planet. The agreement 
is also very good for the 3M,j planet, as long as its eccentric- 
ity is smaller than about 0.1 (see Fig. 8). Planet orbits are 
assumed to always be circular in the one-dimensional models. 



(eo — 0) is much more constant over the course of the sim- 
ulation as its orbital eccentricity remains small (e < 0.02). 
The migration of a 1 Mj with eo = 0.01 proceeds as indicated 
by the short-dashed line in Figure 10. Over the last ~ 1000 
orbits of evolution, however, the migration rate undergoes a 
decrease by a factor ~ 2. During that time, the eccentricity 
grows to 0.09 (see solid line in Fig. 9). 

The effect of non-zero orbital eccentricity on planet migra- 
tion can also be seen in Figure 12, where the evolution of the 
semi-major axis is plotted for simulations with Mp = 1 Afj 
and Mp = 3 A'lj and different initial orbital eccentricities. 
There is a clear trend towards slower migration rates for larger 
orbital eccentricities. In particular, when eo > 0.2, the direc- 
tion of migration is reversed. In all the calculations that show 
outward migration, the angular momentum of the eccentric 
orbit planet increases in time. 

We have conducted some preliminary investigations on the 
cause of this outward torque. One possibility is that it is due 
to the outer disk, as suggested in Papaloizou (2002). When 
the planet eccentricity is large enough, the angular motion of 
the planet at apocenter can be slower than that of the inner 
parts of the outer disk, resulting in form of dynamical friction 
that increases the angular momentum of the planet. For the 
situation we wish to consider, it is not clear how the outer disk 
inner edge is maintained, when this model is applied. This 
disk material loses angular momentum from viscous torques 



and gains angular momentum from the planet, in the usual 
torque balance for a gap. The latter implies that the planet 
should lose, rather than gain, angular momentum. 

A preliminary analysis suggests that the outward torque 
may arise in the coorbital region. This region is supplied 
by material that flows from the outer disk across the gap, as 
discussed in Section 6. In any case, further analysis is required 
to understand this situation. 

We conducted a convergence test on the model with a 1 Mj 
planet and eo — 0.2, which exhibits outward migration. The 
test involved a comparison of the migration rates obtained 
from the grid system GS2 (see Table 2) to those obtained 
from a grid system whose linear resolution was a factor 1.3 
larger everywhere (in both directions on each grid level). How- 
ever, since computing resources were only available to run the 
higher resolution simulation for about 700 orbits, we used the 
Gauss perturbation equations (e.g., Beutler 2005) to compute 
d resulting from the disk's gravitational forces, while keeping 
the planet's orbit fixed. The result of the test is that the 
migration rates, averaged from 200 to 700 orbits, differed by 
only 7% at the two resolution levels. As a check on our use 
of the Gauss equation, we also compared the migration rate 
determined from the Gauss equation, averaged over the last 
100 orbits before release, to the initial d after release, evalu- 
ated by integrating the equations of motion of the planet (see 
the left panel of Figure 12). The two rates, both determined 
on grid system GS2, differed by less than 3%. 

5.3. Effects of Viscosity 

As we discussed in Section 4.3.1, the disk eccentricity de- 
creases with viscosity. For a coupled disk-planet system, we 
similarly expect that the planet eccentricity would decrease 
with Q, since the eccentric corotation resonances become 
stronger. The orbital eccentricity evolution of 3 Mj planets 
in disks with different a values is shown in the top-left panel 
of Figure 13. The model with standard viscosity (solid line) is 
the same as that in Figure 8. Over a period of 2200 orbits, the 
orbital eccentricity of the model with ao = 1.2 x 10"'^ (short- 
dashed line) remains small and never exceeds e ~ 0.01. On the 
other hand, the models with lowest viscosities, qq — 1.2 x 10~^ 
(long-dashed line) and 4x 10^"^ exhibit a generally growing ec- 
centricity. The trend towards faster growth for smaller viscosi- 
ties is confirmed by the model with 1 Mj and ao = 1.2 x 10~^, 
as indicated by the long-dashed line in the bottom-left panel 
of Figure 13. 

The orbital evolution of the semi-major axis of 3 Mj plan- 
ets for the three disk viscosities is displayed in the top-right 
panel of Figure 13. The radial inward migration is faster for 
larger a, as expected in a Type Il-like regime. The plot sup- 
ports the contention that planet eccentricity slows migration. 
Around 4000 orbits after the release time, the two calculations 
with smallest viscosities (qq = 4 x 10~^ and 1.2 x 10"'^) pro- 
duce migration rates respectively equal to d fa —8 x 10~^ ao 
(e f» 0.11) and a « -3x 10"'^ ao (e « 0.14) per orbit. The first 
rate is a factor 8 smaller, while the second a factor 13 smaller, 
than the initial migration speed (i.e., when e ~ 0). For these 
cases, the average eccentricity growth rate at a time of about 
2500 orbits is e « 2 x lO"'' orbit"\ The migration of a 1 Mj 
planet (Fig. 13, bottom-right panel) also shows that, as the 
orbital eccentricity approaches ~ 0.08, |d| starts to reduce. 
When ao = 1-2 x 10"'^, the average migration rate, over the 
last 1000 orbits, is about a factor 5 smaller than it is during 
the first 1000 orbits of evolution after release. 
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Fig. 12. — Left. Semi-major axis evolution of 1 Mj planets for three values of the orbital eccentricity at release: eo = (short- 
dashed line), eo = 0.1 (long-dashed line), and eo = 0.2 (solid line). Right. Same as the left panel but for 3Mj planets with 
orbital eccentricities at release: eo = (short-dashed line), eo = 0.2 (long-dashed line), and eo = 0.3 (solid line). The release 
time is between 1000 and 1200 orbits. Our standard disk parameters were used, including a = 4 x 10""^. 



6. Pulsed Accretion 

For eccentric orbit binary star systems, the accretion from 
a circumbinary disk onto the stars pulsates over the orbital 
period of the binary (Artymowicz & Lubow 1996; Giinther & 
Kley 2002). This effect is related to the pulsating character 
of the equipotential surfaces of the elliptical restricted three- 
body problem (Todoran 1993). 

We measured the mass accretion rate (which we denote as 
Mp) into the inner portion of the planet's Roche lobe (within 
r^cc ~ 0.3 -Rh), following the prescription described in Sec- 
tion 3.2, as a function of the planet's orbit phase. We refer to 
this rate as the planet accretion rate although the flow is not 
resolved on the scale of the planet's radius and thus the rate 
at which the planet would accrete mass may be modulated 
somewhat differently from Mp. Quantity Mp was determined 
by folding the mass accretion rate over planet orbital phase 
and averaging over 500 orbital periods (from 500 to 1000). 
Figure 14 shows the resulting averaged accretion rate, (Mp), 
for 1 Mj and 3M,j planets, versus the true anomaly (i.e., the 
azimuthal position relative to pericenter) of the planet and as 
a function of the orbital eccentricity. The simulations show 
pulsed accretion in cases of eccentric orbit planets. The am- 
plitude of the variability, and to a lesser extent the phase, of 
Mp depends on the orbital eccentricity. 

6.1. Modulation 

The accretion onto a 1 Mj planet with e = 0.1 has two 
asymmetric peaks, the taller of which is around the apocenter 
position (true anomaly equal to n). The secondary peak is 
about 70% of the primary peak and is located close to the 
pericenter position. For e < 0.2, the modulation of Mp in- 
creases with e. For larger orbital eccentricities, modulation 
decreases. This effect may be a consequence of the gap be- 



coming broader and shallower with increasing e (see Fig. 2). 

A similar phase variability is found for the mass accretion 
onto a 3Mj planet (Figure 14, right panel). Even when the 
planet's orbit is circular, (Mp) smoothly varies between 1.5 x 
10"'* Mj and 4.5 x 10"* Mj per orbit. The phasing in this case 
is arbitrary, since the planet's orbit is circular. This behavior 
is related to the eccentricity of the disk. The mass accretion 
is markedly peaked around the apocenter position when e > 
0.1. The mass accretion modulation is again greatest for e — 
0.2. When the planet eccentricity is between e — 0.3 and 
0.4, the highest accretion rate occurs roughly 0.1 orbits after 
apocenter. This delay may be related to the time required by 
material to be captured once it has been perturbed near the 
apocenter. Due to such a delay, accretion on binaries occurs 
near pericenter (Artymowicz & Lubow 1996). 

The density distribution in the vicinity of an eccentric orbit 
planet varies strongly with its orbital phase. This variation is 
illustrated in the panels of Figure 15, which depict the situa- 
tion at the pericenter (left) and apocenter (right), for a 1 Mj 
(top) and a 3 Mj planet (bottom) on an eccentric orbit with 
e = 0.2. The spiral waves have a regular pattern at pericenter, 
when the planet is orbiting in the low-density gap (or cavity). 
As the planet approaches the apocenter, the outer spiral wake 
penetrates higher density regions, which causes fluid elements 
along the wake to lose angular momentum and flow through 
the gap. There are streams of material that extend inwards 
(at r < a and > 0p) which appear in the right panels of 
Figure 15. 

6.2. Mass Growth Timescale 

The mass accretion rate onto a planet with a fixed circular 
orbit decreases with increasing planet mass when Mp > 1 Mj 
(Lubow et al. 1999). The average accretion rate in the simu- 
lations, at t ~ 150 orbits, of a 2 Mj planet is 0.63 times that 
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Fig. 13. — Evolution of planet eccentricity (left) and semi-major axis (right) for 3Afj (top) and 1 Mj (bottom) planets in disks 
with different values of viscosity parameter qq (a at r = ao), ao = 1-2 x 10~^ (long-dashed line), ao — 4 x 10^^ (solid line), and 
ao — 1.2 X 10"^ (short-dashed line). The release time is tj-i^ = 1000 orbits in the top panels and 1100 orbits bottom panels. 



of a 1 A4j planet. The ratio decreases to 0.44 for a 3 Mj planet 
on a circular orbit and at t = 150 orbits. These ratios agree 
within 10% with the values given by Lubow et al. (1999), who 
used an independent code. 

As the disk eccentricity grows, the mass accreted during an 
orbit increases. At later times t > 500 orbits, when S(i_o) ^ 
0.2, the accretion rate onto a 2 Mj planet is 0.71 times the 
rate onto a 1 A/j planet and for a 3 Mj planet is 0.65 times 
the rate onto a 1 Mj planet. Notice that, for a 3Mj planet, 
this implies a 48% increase over the accretion rate at early 
stages, when the disk is circular. These results suggest that 
the eccentricity driven in the disk by a massive planet can 
augment the mass accretion rate onto the planet and hence 
shorten its growth timescale. 

Mass accretion over an orbit period can also be enhanced 



by the planet's orbital eccentricity. For a 1 Mj planet on a 
fixed orbit, the mass growth timescale (defined here as the 
ratio of A'lp to the average accretion rate between 500 and 
1000 orbits) decreased by about 35% when e is increased from 
to 0.2. The reduction of mass growth timescale from e = 
to e = 0.4 is only 17%, perhaps as a result of the wider gap 
and its smoother outer edge at larger orbital eccentricities (see 
Fig. 2). For a 3Mj planet on a fixed orbit with e = 0.3, the 
mass growth timescale is reduced relative to the e = case by 
27%. While for e = 0.4 there was a 13% reduction. For cases 
with e < 0.3, the growth rate was not substantially different 
than the e = case. 

The mass growth timescale of a 1 AIj and a 3 Afj planet 
on circular orbit was also estimated for different values of 
the disk viscosity. For a 3Afj planet an increase in a by a 
factor of 3 over the standard value (see § 2.2) reduced the 
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Fig. 14. — Mass accretion rate in inner parts of the Roche lobe of a 1 Mj (left) and a 3Mj (right) planet as a function of its 
true anomaly and orbital eccentricity. When the true anomaly equals tt, the planet is located at the apocenter of the orbit. 
Orbital eccentricities are listed in the legend of the left panel. The quantity (Mp) is obtained by sampling Afp every 0.02 orbits 
and averaging the outcome from 500 to 1000 orbits. 



growth timescale by about 60%. For both planetary masses, 
a decrease in q by a factor of 3.3 below the standard value, 
lengthened the growth timescale by a factor 2. 

Planetary accretion rates were determined by means of the 
grid system GS2. Simulations executed with grid systems 
GSl and GS3 resulted in very similar outcomes (within 6%) 
for both the modulation of Mp and its average value. The 
results are not sensitive to the smoothing length, e, although 
this parameter can affect the small-scale structure of the flow 
around the planet. We performed a calculation for an e value 
that was reduced by a 25%, with Mp = 3 Mj and e = 0.3, and 
obtained essentially the same accretion rates (within 1%). 

6.3. Accretion towards the Star 

The region interior to a planet's orbit likely contains an 
inner disk which cannot be resolved in the current two- 
dimensional simulations (see discussion in § 5.2). We esti- 
mate the modulation of mass onto this disk as a function of 
orbital phase by considering the mass flow rate across the in- 
ner boundary. As material accretes through the inner disk, 
the modulation would be expected to weaken. It is unclear 
whether the modulation would be reflected as a variable accre- 
tion at the surface of the star. Perhaps it could be manifested 
as variability in emission from the region where the inflow 
meets the outer edge of the inner disk. 

Figure 16 plots the accretion rate at the inner boundary 
(A/*) versus the true anomaly of the planet for cases with 
A/p — 1 Mj and 3 Afj . The case involving the more massive 
planet produces an accretion rate through the inner bound- 
ary that is largest when the planet is close to the apocenter. 
The maximum of {M,} occurs before the apocenter passage 
for the case involving the 1 Mj planet. The average mass ac- 
creted by the star during one orbital period of the planet is 
5.8 X 10"* A/, and 9.9 x 10~*A//. for the 1 Af j and the 3Afj 
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Fig. 16. — Mass accretion rate towards the star at the in- 
ner boundary r = 0.3 ao as function the true anomaly of the 
planet. The mass accretion rate is averaged over several planet 
orbital periods at t ~ 1000 orbits. The dashed and solid 
lines refer to models with A^p = 1 Afj and 3 A/j , respectively. 
The planet's orbital eccentricity is e = 0.2. When the true 
anomaly is tt, the planet is located at the apocenter. 



cases, respectively. The same models with no orbital eccen- 
tricity show relatively constant (Af*) values of 1 x 10~^Af» 
(A//p = lAfj) and 6.3 x 10"* (Afp = 3A^j) per orbit. As a 



G. D'Angelo, S. Lubow, & M. Bate 



17 



-&.oa 



-5.B1 



-5JZ2. -A.B3 



-i.44 



-&.oa 



-4.50 



-3.75 



-3.no 



3,30 -: 



3.20 - 



3,10 -i 



3,00 T 



2.30 




3,40 



3,30 



3.20 



3,10 



3,00 



2. SO 




'I II III II 1 1 III II II 1 1 1 II II III II 1 1 II III II 1 1 III II II 1 1 1 II II III II' 

1.10 1.20 1.3D 1.40 



-h.Oa -5.BD -5JZ0 -4.79 -4.39 



-h.Oa -5.B1 -5JZ2, -4.B4 -4.4S 



3,40 




2. so -' 



3,30 -: 




2. so - 



0.6 0.7 D.a 0,9 



1,00 1.10 1.20 1,30 1.40 



Fig. 15. — Density structure around a 1 Af j (top) and a 3Afj (bottom) planet at pericenter (left) and apocenter (right). The 
vertical axis is the azimuth about the star and the horizontal axis is the distance from the star in units of ao. The orbital 
eccentricity is e = 0.2 and t ~ 1000 orbits. The solid circle indicates the instantaneous location of the planet. A surface density 
of 10^^ corresponds to 3.29 gcm^^. 



comparison, the mass accretion through a steady-state a-disk 
is SttE u (Lynden-Bell & Pringle 1974; Pringle 1981) which, 
using the initial unperturbed surface density and standard 
viscosity (qq = 4 x 10~^), yields 2.5 x 10"'' M« per orbit or 
2 X 10"** A'fgyr^V 

The ratio of the accretion on the star to the accretion in 
the disk, outside the gap, can be expressed as (Aft)o/((Af*)o + 
(Afp)o), where the subscript "o" denotes is the integral of the 
respective accretion rate over a planetary orbit. For both 
the 1 Mj and the 3Afj planets on circular orbit, this ratio is 
0.19. When e = 0.2, the ratio is 0.09 and 0.26 for the 1 j 
and the 3 Mj planet, respectively. The reduced mass transfer 
across the planet's orbit, when Mp = 1 Afj and e — 0.2, can 
be attributed to the increased accretion rate onto the planet. 



as reported above. In the 3 Mj case, {Mp)o does not vary 
significantly as e varies from to 0.2. Instead, the mass flux 
across the gap is likely enhanced by the radial excursion of 
the planet (see § 6). 



7. Summary and Discussion 

We simulated the orbital evolution of circular and eccen- 
tric orbit giant planets embedded in circumstellar disks. The 
disks were analyzed using a two-dimensional hydrodynamics 
code that utilizes nested grids to achieve high resolution in a 
large region (2 a x 27r/3) around the planet. The disks were 
modeled as an a-disk and a few values of a were considered. 
We investigated planet masses of 1 Mj , 2 Afj , and 3 Afj and 
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initial orbital eccentricities that ranged from to 0.4. 

Disk gaps become broader and shallower as the planet ec- 
centricity increases (see Fig. 1 and 2). The density near the 
orbit of the planet is very small compared with the density in 
the disk for all eccentricities considered. A planet on a fixed 
circular orbit can cause an initially circular disk to become 
eccentric (see Fig. 4). The disk eccentricity is suppressed at 
lower planet masses (Mp < 1 Mj) and higher disk viscosities 
(a > 0.01), as also found by Kley & Dirksen (2006) and by 
Papaloizou et al. (2001) at higher planet masses. We attribute 
the eccentricity growth to a tidal instability associated with 
a series of eccentric outer Lindblad resonances in the inner 
parts of the outer disk (Fig. 7). The same type of instability, 
involving an inner disk, is thought to be responsible for the 
superhump phenomena in binary star systems (Lubow 1991a; 
Osaki 2003). 

The simulations indicate that planet eccentricity can grow, 
as a consequence of disk-planet interactions (Fig. 8 and 9). 
The growth is stronger in the 2 Mj and 3 Mj cases than for 
IMj, and for lower disk viscosity (q < 4 x 10"^). Planet 
eccentricities of --^ 0.1 were found in the simulations over the 
course of a few thousand orbits for 2 Mj and 3 Afj planets. 
A similar eccentricity growth is obtained for a 1 Mj planet in 
a disk with viscosity a ~ 10""^. The planet and disk both 
acquire eccentricity as they interact, which may lead to com- 
plicated time-dependent behavior of their eccentricities. The 
planet eccentricity growth is likely aided by the disk eccen- 
tricity growth. The results suggest that the eccentric growth 
found for 10 Mj planets by Papaloizou et al. (2001) also oc- 
curs for lower planet masses. The higher resolution achieved 
by our calculations may be playing a role in obtaining this 
growth. 

For circular orbit planets, migration occurs on roughly the 
local viscous timescale, as expected for Type II migration. 
However, it is slowed for eccentric orbit planets. This result 
appears for several configurations with either dynamically de- 
termined (Fig. 10 and 13) or imposed planet eccentricities 
(Fig. 12). For a 2Mj case, even migration reversal (outward 
migration) is found for a dynamically determined eccentricity 
(Fig. 10). Migration slowing or reversal would have important 
consequences for the planet formation process. The cause is 
not yet clear. It may involve torques from outer disk (Pa- 
paloizou 2002) or instead from the coorbital region. Some 
preliminary evidence suggests the latter. 

Mass accretion both within a planet's Roche lobe and 
through a gap can be strongly modulated with orbital phase 
for eccentric orbit planets or eccentric disks (Fig. 14, 15, 
and 16). The modulation was largest for planet eccentricity 
e ~ 0.2. This pulsating accretion is similar to what is found 
for eccentric orbit binary stars embedded in a circumbinary 
disk (Artymowicz & Lubow 1996), although the phasing is 
different. Both disk and planet eccentricity also lead to en- 
hanced accretion onto the planet. This enhancement likely 
helps planets achieve higher masses. 

The simulations lend support to the idea that disk-planet 
interactions cause planet eccentricity growth, along the lines 
of Goldreich & Sari (2003). The simulations suggest that 
planet eccentricities are easier to achieve for higher mass plan- 
ets (Mp > 2Mj). Our results are subject to the usual limita- 
tions in approximate initial conditions, simulation time, radial 
range for coverage of the disk (likely resulting in the lack of an 
inner disk), the a-disk model, and the use of various numer- 
ical devices. We also neglected disk self-gravity, which may 



affect migration especially for higher mass disks (Nelson & 
Benz 2003a). 

However, it is not clear that typical extra-solar planet ec- 
centricities of 0.2-0.3 can be achieved through disk-planet in- 
teractions. The eccentricity growth at later times shows in- 
dications of slowing and possibly stalling for e < 0.15 (see 
Fig 13.) Perhaps higher eccentricities can be achieved for disks 
with different properties (e.g., lower viscosity and smaller 
disk's aspect ratio). Eccentricity may be limited by damping 
due to high order eccentric inner Lindblad resonances that lie 
outside a planet's orbit. Simulations of eccentric orbit binary 
star systems suggest that little eccentricity growth occurs for 
e > 0.5 (Lubow & Artymowicz 1992). Although the simu- 
lated planets do not achieve orbital eccentricities in excess of 
0.15 over the duration of the simulated evolution (for config- 
urations that start from circular orbits), the simulation times 
correspond to less than 10* years. Migration slowing and re- 
versal may permit the planets to achieve higher eccentricities 
on longer timescales while avoiding orbit decay into the disk 
center/host star. 

We have not yet investigated eccentricity evolution of sub- 
Jupiter mass planets. They may also provide important con- 
straints. Other simulations suggest that disks with standard 
viscosity have only mild gaps for smaller planet masses of 
Mp < 0.1 Mj (e.g., D'Angelo et al. 2003; Bate et al. 2003). 
Under those conditions, disk-planet interactions likely lead to 
eccentricity damping, due to the dominance of the coorbital 
Lindblad resonance (Ward 1986; Artymowicz 1993). The ob- 
servational determination of eccentricities for small mass plan- 
ets would help constrain these models. The planet around 
HD 49674 is close to this regime. It has a minimum mass of 
0.11 Mj and a best-fit eccentricity of 0.29 (P. Butler, private 
communication). Since it is close to the central star (the pe- 
riod is 4.9 days), it is possible that the eccentricity evolution 
could be more complicated, especially if it became trapped in 
a central disk hole. Examples of isolated planets like this, but 
at longer periods would provide useful constraints. 
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